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ABSTRACT 

We present a sub-50 pc-scale analysis of the gravitational lens system SDR81 at redshift 
3.042 using Atacama Large Millimetre/submillimetre Array (ALMA) science verification 
data. We model both the mass distribution of the gravitational leasing galaxy and the pix¬ 
elated surface brightness distribution of the background source using a novel Bayesian tech¬ 
nique that fits the data directly in visibility space. We find the 1 and 1.3 mm dust emission to 
be magnified by a factor of //tot = 12.6 ± 0.4, giving an intrinsic total star-formation rate of 
315 ± 60 Mq yr“^ and a dust mass of 6.4 ± 1.5 x 10^ Mq. The reconstructed dust emission is 
found to be non-uniform, but composed of multiple regions that are heated by both diffuse and 
strongly clumped star-formation. The highest surface brightness region is a ^1.9 x 0.7 kpc 
disk-like structure, whose small extent is consistent with a potential size-bias in gravitation¬ 
ally lensed starbursts. Although surrounded by extended star formation, with a density of 
20-30 ±10 Mq yr“^ kpc“^, the disk contains three compact regions with densities that peak 
between 120-190 ±20 Mq yr“^ kpc“^. Such star-formation rate densities are below what 
is expected for Eddington-limited star-formation by a radiation pressure supported starburst. 
There is also a tentative variation in the spectral slope of the different star-forming regions, 
which is likely due to a change in the dust temperature and/or opacity across the source. 
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1 INTRODUCTION 

Submillimetre (sub-mm) galaxies are understood to be the main 
drivers of star-formation activity at high redshift and are believed 
to be progenitors of present-day early-type galaxies (e.g. |Blain et al.| 
|2002[ |S winbank et al.|2006| l. Their large far-infrared luminosities 
10^^ Lq) are due to stellar ultraviolet (UV) radiation being ab¬ 
sorbed and re-emitted at far-infrared wavelengths within dusty star¬ 
forming regions. Thus, studies of the dust emission at (rest-frame) 
far-infrared wavelengths can be used to trace star-formation activity 
directly, even if the stars themselves suffer from heavy extinction 
due to the dust ( |Kennicutt|1998| l. 

While detailed studies of the star-formation processes within 
these objects have not been possible due to the arcsecond-scale an¬ 
gular resolution of the early instruments (e.g. SCUBA; MAMBO), 
the advent of interferometric arrays, such as the Atacama Large 
Millimetre/submillimetre Array (ALMA) and the Submillimetre 
Array (SMA), has led to a significant improvement in our knowl¬ 
edge about the structure of starbursting galaxies at these wave¬ 
lengths ( [Karim et al.|2013| ). Furthermore, the study of gravitation- 
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ally lensed sub-mm galaxies has allowed an investigation of the 
formation of high redshift galaxies at a resolution comparable to 
what has already been achieved at lower redshift ( |S winbank et al.j 
|2010j ). Due to the combination of the magnification provided by the 
foreground gravitational lenses, and the fact that the observed flux 
density of sub-mm sources is roughly constant across a wide range 
of redshifts (the so-called negative K-correction), distant sources 
up to a redshift of 5.6 have been detected in recent sub-mm 
surveys (jNegrello et al.|2010l|Bussmann et al.|2012|| Vieira et al.j 

|20T^ . 

It is now possible to study the detailed star-formation pro¬ 
cesses within these galaxies at <100 mas angular resolution with 
the long baselines provided by ALMA. Such studies will allow the 
size and structure of the heated dust emission to be directly ob¬ 
served for the first time, which is important for determining the 
actual magnification of the emission regions without the effects 
of differential magnification. This can establish whether there is 
a size-bias in gravitationally lensed starbursting galaxies and place 
robust limits on the star-formation rate density within these galax¬ 
ies. 

In this letter, we present the analysis of high angular res¬ 
olution imaging of the heated dust continuum emission from 
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Figure 1. The ALMA continuum imaging of SDR 81 at 236 GHz (Upper) 
and 290 GHz (Lower). The images have been produced at two angular reso¬ 
lutions using a cut in the wv-data of 0.5 (left) and 2 MA (right), and weight¬ 
ings (natural and Briggs robust = 0) to emphasise the structure seen on dif¬ 
ferent scales. The synthesised restoring beam is shown as the white ellipse 
in the bottom left corner of each map. These maps are only shown to illus¬ 
trate the surface brightness distribution of the source. It is the visibility data 
that are used for the gravitational lens modelling and source reconstruction. 


the gravitationally lensed starburst galaxy SDP.81 (H-ATLAS 
J090311.6+003906; g = 3.042± 0.001; [Negrello et al.||20T0| 


Prayer et al.|201 l||Bussmann et al.|2013||Dye et al.|2014||Negrello| 
et al.|2014> that were taken during the science verification phase of 
ALMA. These observations were taken using the new long baseline 
capability of ALMA and provide a unique view of this star-forming 
galaxy. We determine the structure of the background source by ap¬ 
plying a Bayesian pixellated visibility-fitting lens modelling tech¬ 
nique ( [Rybak, Vegetti & McKean||2015| ). A reconstruction of the 
emission line properties of the background source, using the lens 
model derived here is presented in a companion paper (Paper II). 


2 ALMA CONTINUUM IMAGING 

The ALMA science verification dataset of SDP.81 was taken over 
several observing blocks during 2014 October and November. The 
data were taken in Bands 4, 6 and 7 and comprised both continuum 
imaging datasets centred at 151,236 and 290 GHz, and spectral line 
datasets centred on the CO (5-4), CO (8-7), CO (10-9) and H 2 O (2- 
1) emission lines. In this letter, we only consider the datasets taken 
in Bands 6 and 7 since they have the highest continuum signal-to- 
noise ratio. 

The Band 6 and 7 continuum datasets were each taken using 
2 spectral windows of 2 GHz total-bandwidth and with 128 chan¬ 
nels in each of the linear polarisations (XX and YY). Other spectral 
windows were available, but these were used primarily for spectral 
line studies. A visibility averaging time of either 2 or 6 s was used. 
The observations were carried out with 31-36 antennas in the array 
and baseline lengths from ~ 15 m to 15 km, with 10 per cent of the 
baselines shorter than 200 m. This gave an array that was sensitive 
to structures on angular scales between 19 arcsec and 16 mas (at 


236 GHz). We refer to the ALMA partnership (2015) for more de¬ 
tails about the observations and data processing. We found that the 
visibility data for the final observing block of the Band 6 dataset 
had a larger rms noise, and so this observing block was removed 
from the dataset. Also, the absolute fiux-density calibration of the 
first spectral window of the Band 6 dataset was higher by ~28 per¬ 
cent relative to the other sub-bands, an off-set that we corrected for 
during our analysis. The post-reduction data quality of the Band 7 
dataset was excellent. The total time on-source was 4.4 and 5.6 h 
for the Band 6 and 7 observations, respectively. 

The main goal of this letter is to reconstruct the surface bright¬ 
ness distribution of the background source from the visibility data. 
However, each dataset had over 10® visibilities, which is computa¬ 
tionally expensive to analyse. To increase the speed of our analysis, 
we therefore averaged each spectral window to a single channel and 
used a visibility averaging time of 20 s, both of which have minimal 
bandwidth and time smearing effects on the final dataset. In addi¬ 
tion, from our imaging tests, we found that the recovered structure 
of the source is highly dependent on the baseline data that are used. 
In Fig.we present the imaging data for Bands 6 and 7 using a cut 
in the wv-data at baseline lengths of ^ 0.5 and ^ 2 MA (where A is 
the observing wavelength; 1.3 and 1 mm at Band 6 and 7, respec¬ 
tively). The ^0.5 MA images show that at mm-wavelengths the ob¬ 
served system comprises four images in a cusp configuration with 
evidence for a low-surface brightness Einstein ring. The ^ 2 MA 
image shows that the large arc is composed of several structures 
with varying surface brightness, and although the Einstein ring is 
effectively resolved out at these scales, there is still evidence of 
emission from this part of the source. The total fiux-density of the 
source at Bands 6 and 7 at ^ 0.5 MA from our fit to the visibility 
data are 21.7 zb 2.2 and 38.9 zb 3.9 mJy, respectively, assuming a 10 
percent uncertainty in the fiux-density calibration. 

Eor the wv-data at ^ 1 MA (i.e. removing the shortest base¬ 
lines), we found that the extended arc is mostly resolved out, which 
demonstrates that most of the emission at mm-wavelengths is ex¬ 
tended; it is only the most compact components within the indi¬ 
vidual images that are detected with the > 2.5 km baselines. We 
carried out additional tests to determine the scales on which most 
compact structures are resolved out. We found that they are de¬ 
tected up to around 5 MA, that is, around 30 mas-scales. Based on 
these results, we used the visibility data at wv-distances ^ 2 MA 
for our lens modelling only, because there is little or no extended 
structure detected on baselines longer than this. Also, a cut was 
used instead of a wv-taper because this limited the number of visi¬ 
bilities needed for the gravitational lens modelling to a manageable 
dataset of about 1.1-1.5 x 10^ visibilities per spectral window. 

With a Briggs weighting scheme (robust = 0), which best 
represents the intrinsic weighting scheme of visibility lens mod¬ 
elling within a Bayesian framework ( jJunklewitz et al.|2013[|Rybak| 
jet al.|2015] l, the synthesised beam-size of the ^ 2 MA dataset was 
95 X 71 mas at a position angle of 64.5 deg east of north (at 236 
GHz). We note that the emission line data are also only detected 
on the shorter baselines, and so our choice of wv-cut will allow the 
continuum and molecular gas emission to be compared on similar 
angular-scales (see Paper II). 


3 GRAVITATIONAL LENS MASS MODELLING AND 
SOURCE RECONSTRUCTION 

In order to reconstruct the smooth lensing potential and the back¬ 
ground source surface brightness distribution, we use the Bayesian 
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Figure 2. (Left) The image-plane residuals from the Fourier transform of the residual visibilities (data—model) and gridded assuming uniform weighting, 
(left-middle) the modelled image-plane brightness distribution with critical curves, (right-middle) the reconstructed source-plane brightness distribution with 
caustics, (right) and the uncertainty on the reconstructed source surface brightness distribution at Band 6 (upper) and at Band 7 (lower). The image-plane 
residuals are in units of mJy beam“^. The peak surface brightness of the reconstructed source is normalised to 1 units pixel” ^. The lens centre is at position 
(0, 0). The source surface brightness crosses the caustics, which results in parts being either quadruply or doubly imaged. There is also diffuse and extended 
structure that produces the Einstein ring emission. 


visibility-fitting method of [Rybak, Vegetti & McKean| ( |2015| ). This 
method extends the gravitational imaging technique introduced by 
[Vegetti & Koopman^ ( |2009j ) to interferometric data. 

Due to the incomplete sampling of the wv-plane, the recon¬ 
struction of the sky surface brightness distribution from the mea¬ 
sured visibility function is a non-trivial inverse problem. More¬ 
over, the pixel-to-pixel noise in the sky surface brightness maps 
is correlated, and since the individual baselines filter structure in 
the source on various angular-scales, the observed surface bright¬ 
ness (in coarsely resolved imaging) is not necessarily conserved. 
As a result, fitting to the data directly in visibility space has be¬ 
come the established method to determine the lens model param¬ 
eters and to reconstruct the surface brightness distribution of the 
background source from interferometric data (e.g. [Wucknitz|2004| 
[Hezaveh et aL|201^[Bussmann et al.|2013[ l. However, these anal¬ 
yses have so far been dependent on assuming a parametric form, 
typically multiple Gaussian or Sersic profiles, for the background 
source. These are unlikely to represent well the complex surface 
brightness distribution of dusty star-forming galaxies on sub-kpc 
scales ( [Rybak et al.[2015[ ). It is for this reason that we instead use a 
pixelated source model to reconstruct the extended emission from 
the high angular resolution ALMA observations. 

The technical implementation of pixelated source grids and 
their application to interferometric data are discussed in detail by 
[Vegetti & Koopman^ ( [2009[ ) and [Rybak et al.[ ( [20T5] t, respectively. 
To summarise, the source-plane grid is constructed using a De¬ 
launay tessellation so that the highly magnified regions are sam¬ 
pled more densely than those with lower magnification, thus sav¬ 
ing computational time. This also enables the small-scale structure 
in the highly magnified regions of the source to be better sampled. 
The modelling method optimises the Bayesian evidence by min¬ 
imising a combination of the measured in the visibility space 
and a regularisation performed on the source plane grid. By includ¬ 
ing the regularisation term, we break the ill-posed nature of the 
lensing equation by enforcing a degree of smoothness in the source 
surface brightness distribution, while avoiding overfitting the noise. 
Here, we adopt the gradient rather than the curvature regularisation 


to prevent the source from being overly smoothed. A detailed dis¬ 
cussion of various levels of source regularisation is given by[Suyu[ 

Fnri ( [20061 ). 

We use an elliptical power-law mass distribution with an exter¬ 
nal shear to describe the lensing potential. This model is defined by 
a mass parameter (/^o; the normalization of the surface-mass den¬ 
sity), an ellipticity (q) and position angle (0; east of north), a power 
law density-profile (7), a shear (L) and a shear position angle (L^; 
north of west). The most probable a posteriori model parameters 
are presented in Table with their uncertainties. These are con¬ 
sistent with the gravitational lens modelling of the SMA data by 
[Bussmann et al.[ ( [2013[ ) and of the HST data by [Dye et aL\ ( |2014[ ) 
at the 2cr-level, with the small differences being attributed to the 
choice of the mass model, for example, Bussmann et al. do not in¬ 
clude shear in their model. 

With the lens model in hand, we reconstruct the source sur¬ 
face brightness map for each of the two spectral window datasets 
in Bands 6 and 7 separately, and then take an average to improve 
the signal-to-noise ratio. This is done partly to provide a consis¬ 
tency check between the source reconstructions for each spectral 
window, but also to limit the size of the visibility datasets that are 
fitted. We estimate the noise level for each baseline at a given time 
and frequency by calculating the rms of the visibilities measured 
within the same observing block (typical length of approximately 
30 minutes). This could lead to a biased noise estimate if the sys¬ 
tem temperature varies sharply on short time-scales, but this was 
not found to be the case. 

In Fig. [^ we show the average reconstructed source surface 
brightness distribution and the resulting sky-model for the Band 6 
and 7 datasets. In addition, we also present error estimates for the 
reconstructed source. These are taken from the rms scatter in the 
source reconstructions between the different spectral windows of 
each band. We find that the structure of the source surface bright¬ 
ness distribution is in general consistent in the two bands, which 
provides an independent check of our source reconstruction. How¬ 
ever, there are variations in the relative surface brightness of the 
structure within each band (see below for discussion). Note that the 
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high surface brightness inner region of the source is reconstructed 
at higher resolution and at a higher signal-to-noise ratio compared 
to the rest of the source plane. 

We have determined the residuals between the data and the 
model in the visibility plane, and find that they are within the un¬ 
certainties of the individual visibilities in the real and imaginary 
parts. The residuals are also within the expected rms over the whole 
dataset. For example, 4.7 and 0.3 percent of the visibilities are out 
with the 2- and 3-cr levels, respectively, which is in good agree¬ 
ment with what is expected from a Gaussian distribution (white 
noise). However, to illustrate the goodness of fit of our model, we 
also show in Fig.j^the image plane residuals. Note that the image 
plane residuals are dependent on the weighting scheme used when 
gridding the visibilities. 


4 THE INTRINSIC SOURCE PROPERTIES OF SDP.81 

The reconstructed continuum emission of SDP.81 at 236 and 290 
GHz possesses a complex morphology; there are three distinct re¬ 
gions of the source that are resolved in our pixelated reconstruc¬ 
tions, but there is also evidence of a wider region of low-surface 
brightness emission (see Fig. |^. This outer region is consistent 
with the extended emission that forms the Einstein ring seen in 
our low resolution images (see Fig.[^. It is clear that the sub-kpc 
structure of the galaxy is not uniform, but contains several regions 
of intense dust emission that are presumably associated with re¬ 
cent star-formation. In particular, the central part of the source is 
elongated and has an extent of ~1.9x0.7 kpc. We interpret this 
as evidence of a star-forming disk (e.g. as in GN20; [Hodge et al.| 
|2015j ), but it could also be due to the multiple components of an 
ongoing merger (e.g. jSwinbank et al.||2006| [Tacconi et al.||2008l 
[Engel et al.|2010| |. We note that photo-dissociation region (PDR) 
models for SDP.81 from an analysis of the [O III] lines in the far- 
infrared spectrum, as measured with Herschel, constrain the size of 
the emitting region to be 500-700 pc (effective radius; [Valtchanov[ 
[et al.|20lT[ l, which is consistent with the overall size of the central 
star forming region that we detect. However, the individual regions 
with the most intense dust emission are smaller than this size, that 
is, ^ 500 pc. 

The extent and surface brightness distribution of the recon¬ 
structed source are considerably different than reported by [Buss-[ 
[mann et ^ < |2013[ ), who modelled the SMA data at 880 jim (340 
GHz) using a single Sersic profile with an effective radius of 
~4.1 kpc and at a position angle that is almost perpendicular to 
the extension seen in the ALMA reconstruction. However, as [Buss-[ 
[mann et si] ( [2013[ ) discuss, their single parametric source model 
is too simple to fit the observed surface brightness distribution of 
the 880 /xm emission at 0.6 arcsec angular resolution; this high¬ 
lights the advantage of using pixelated source reconstructions and 
100 mas angular resolution observations to study the intrinsic prop¬ 
erties of gravitationally lensed starburst galaxies. 

We estimate the lensing magnification /x of the source 
by calculating the ratio of the reconstructed emission in the 
source-plane and image-plane. The total magnification is found 
to be /xtot = 17.6 zb 0.4. This magnification is dominated by 
the central star-forming disk, which has a magnification of 
MSF disk = 25.2 zb 2.6. Our total magnification is a factor of 1.6 
higher than those reported by [Bussmann et al.[ ( [2013[ ) and [Dye[ 
[et al.[ ( [20T4] ), which is most likely due to the differing source struc¬ 
ture that we find compared to their studies; compact sources that 
are located close to the lens caustic can have very high magnifi¬ 


cations. In addition, our analysis directly measures the magnifica¬ 
tion of the resolved continuum emission from the dust, which is 
most relevant when studying the star-forming properties of SDP.81. 
The compact structure of the source and the high magnification that 
we find are consistent with a potential size-bias in gravitationally 
lensed sources (e.g. [Tacconi et al.[[2008] l and the results for this 
source should not necessarily be extrapolated to the general sub- 
mm galaxy population. 

Using our measurement of the total magnification, we have 
determined the intrinsic properties of SDP.81 from a fit to 
the UV/optical-to-sub-mm broad-band spectrum that was car¬ 
ried out by [Negrello et al.[ ( 2014| ). We find a star-formation 
rate of 315 zb 60 M© yr“\ a far-infrared luminosity of 
3.1 zb 0.4 X 10^^ L© and a dust mass of 6.4 zb 1.5 X 10® M©. Here, 
we only consider the properties derived from the far-infrared part 
of the spectrum, since this is the parameter space where our magni¬ 
fication is valid. These results are consistent with what is typically 
found for sub-mm galaxies at redshift ~2.5 (e.g. ~400 M© yr“^; 
[Coppin et al.|2008[ ). 

In Pig. we show the star-formation rate density of the re¬ 
constructed source, under the assumption that all the heated dust 
emission is solely due to star formation. Also shown is the un¬ 
certainty in the star-formation rate density, which was calculated 
from the rms over all spectral channels. We find that the ex¬ 
tended component has a moderate star-formation rate density of 
~20-30 M© yr“^ kpc“^, but there are multiple regions of in¬ 
tense star-formation (^100 M© yr“^ kpc“^) as opposed to a sin¬ 
gle starburst site. The star-formation disk is forming stars at a 
mean rate of over 100 M© yr“^ kpc“^. This is several orders 
of magnitude higher than is seen within disk galaxies in the lo¬ 
cal Universe ( [Leroy et al.[[2013] ), but is comparable to the ex¬ 
tended disk seen in GN20 at redshift 4 ( [Hodge et al.[[2015| l. The 
maximum star-formation rate density in our reconstructed map is 
190 zb 20 M© yr“^ kpc“^, which is below the theoretical expec¬ 
tation for Eddington-limited star-formation by a radiation-pressure 
supported starburst (~1000 M© yr“^ kpc“^). There are several 
other examples of high redshift sub-mm luminous galaxies that 
have sub-Eddington starbursts ( [Bussmann et al.|2()T^ [Hodge et al.[ 
[2015| ), but this is the first time that the star-formation density of a 
galaxy has been directly mapped on sub-50 parsec-scales thanks to 
gravitational lensing and the new high resolution imaging capabil¬ 
ity of ALMA. Clearly, high resolution imaging of the larger sam¬ 
ple of gravitational lenses found during the SPT and H-ATLAS lens 
surveys, coupled with pixelated source reconstructions, is needed to 
determine if these galaxies as a population are considerably below 
the Eddington-limit for radiation-pressure supported starbursts. 

We also see tentative evidence for a change in the dust con¬ 
tinuum slope over the resolved extent of the intense star-forming 
regions within SDP.81 (see Pig. for the relative brightnesses of 
the regions). Between 230 and 290 GHz, the flux-ratio of the en¬ 
tire source is 0.56 zb 0.08 (or a spectral index of a = 2.8, where 
Su oc zx^), which is as expected for the best fit modified black 
body model for the dust emission (temperature Td = 34 zb 1 K and 
emissivity [3 = 1.5; [Bussmann et al.|2012| l. We find that the cen¬ 
tral dominant clump of dust emission in our maps has a flux-ratio 
of 0.51 zb 0.07 between 230 and 290 GHz, but the two additional 
clumps of intense dust emission to the north and south have fiux- 
ratios of 0.41 dz0.03 and 0.45 zb 0.05, respectively (equivalent to 
a ~ 4.3-3.9). This change in the spectral slope of the continuum 
emission suggests that the different sites of intense star-formation 
have different physical conditions, but further observations at high 
angular resolution with Band 9 (and 10) will be needed to deter- 
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Table 1. The best lens model parameters for SDP.81 derived from an analysis of the ALMA continuum visibility data. 


KQ (arcsec) 

q 

0 [deg] 
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r 

re [deg] 

L606±0.005 

0.82±0.01 

8.3±0.4 

2.00±0.03 

0.036±0.004 

3.0±0.2 




A X [kpc] 


Figure 3. (Upper) The star-formation rate density (left) and the un¬ 
certainty (right) of the SDP.81 reconstructed source at redshift 3 
shows multiple regions of intense star-formation. The colour-bar is in 
units of Mq yr“^ kpc“^ and the peak star-formation rate density is 
190 ±20 Mq yr“^ kpc“^. (Lower) The flux-ratio between 236 and 290 
GHz (left) and the uncertainty (right) for a region with a star-formation rate 
density ^ 80 Mq yr“^ kpc“^. The colour-bar is in units of flux-ratio. For 
all plots the iso-contours of the star-formation rate density are set at 20, 40, 
80, 120 and 160 Mq yr-^ kpc-^. 

mine the extent of the variation in the dust temperature and/or opti¬ 
cal depth at these locations. However, our results show that global 
temperature models are likely too simple to explain the emission 
from galaxies with multiple sites of star-formation, as expected. 


5 SUMMARY 

We have presented new continuum imaging of the gravitationally 
lensed starburst galaxy SDR81 at redshift 3.042 using the new 
long baseline capability of ALMA. We find that on ~0.5 arcsec- 
scales, the background source at mm-wavelengths is gravitationally 
lensed into four images in a standard cusp configuration with a low- 
surface brightness Einstein ring. On ~100 mas-scales, the back¬ 
ground source is resolved into two gravitational arcs that clearly 
show multiple regions of enhanced dust emission. Most of the 
source is completely resolved out on 30-mas-scales, which corre¬ 
sponds to baseline lengths of MA. Our pixelated source recon¬ 
struction from the visibility data, shows that the unlensed source is 
~1.9 kpc in diameter and elongated; consistent with a star-forming 
disk that has multiple regions of intense star-formation. This small 
linear size of the source is consistent with the expected size-bias as¬ 
sociated with gravitationally lensed starbursting galaxies. The star- 
formation rate density varies from ~20 to 190 Mq yr“^ kpc“^ 
across the source, which is as expected for a sub-Eddington-limited 
starburst. We also detect a tentative variation in the far-infrared con¬ 
tinuum slope across the source. 


The results presented in this letter represent the first direct 
mapping of the star-formation density of a galaxy at such small 
scales and highlight the power of combining gravitational lensing 
and the high resolution imaging capability of ALMA. 
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